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1. Introduction 

The creation and motion of large numbers of crystal lattice dislocations is the most 
fundamental feature of crystal plasticity. During the last half century, the physical 
properties of individual dislocations and their interactions with localised obstacles have 
been studied extensively On the other hand, the complex collective dynamics of strongly 
interacting many-dislocation systems is still far from being understood. Fortunately, 
everyday plastic deformation processes very often proceed orders of magnitude slower 
than the typical relaxation times of the underlying dislocation system. These conditions 
often permit to study the problem in a quasistatic approximation [TJ [2]. Beyond the 
quasistatic limit, however, much less work has been devoted to studying the dynamics of 
collective dislocation motions which lead to the formation of metastable configurations, 
and to transitions between such configurations in driven dislocation systems. However, 
such collective motions are crucial for understanding rapid dislocation processes which 
not only occur in shock loading but, in the form of dislocation avalanches, are a generic 
feature of the dynamics of driven dislocation systems [3]. 

The first studies of dynamic relaxation processes in dislocation systems were 
performed by Miguel et al. with the protocol of applying a constant external shear stress 
to well relaxed dislocation configurations [U [5j [6]. The ensuing creep relaxation was 
numerically shown to follow Andrade's law stemming from the underlying intermittent 
and correlated motion of dislocation structures. The connection between the mesoscopic 
and macroscopic features of the process was, however, not analysed in detail. 

Another direction was taken by the present authors who conducted systematic 
studies of the relaxation dynamics of initially random configurations of straight 
dislocations. This is an important issue since the elastic energy density £ of a random 
dislocation system of density p is known to diverge with the logarithm of system size 
L, £ oc pln(L/b) [7\ where b is the modulus of the dislocation Burgers vector. In 
a well-relaxed dislocation arrangement, on the other hand, the same quantity scales 
like £ oc pln(l/(6 v /p)), i.e., the screening length corresponds to the mean dislocation 
spacing [8j [9j [10]. As the mean square stress is proportional to the elastic energy 
density, this screening also removes a logarithmic divergence of the width of the internal 
stress probability distribution [TT], and of the X-ray line width [101 IH EG] . Numerical 
experience showed that, at least in single slip geometries, the relaxation processes that 
lead to screened dislocation arrangements exhibit slow, power law characteristics for 
quantities such as the elastic energy or the average dislocation velocity [13]. A model was 
proposed which relates the power-law relaxation dynamics to the gradual extinction of 
initial dislocation density fluctuations [14J . The present paper presents a comprehensive 
numerical investigation which allows to check in detail the model predictions and 
complements the earlier work by extending the investigation to multiple slip geometries 
and to dislocation systems of non-zero net Burgers vector, and by studying the influence 
of an external driving stress on the relaxation process. 

The paper is organised as follows. In section [2] the problem is defined and technical 
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details of the simulations are presented. Section [3] unfolds a scaling model of the 
relaxation process from a chemical analogy and uses this model to predict the evolution 
of simulation measurables. Section H] then gives a detailed comparison between model 
predictions and numerical results. The results are discussed and conclusions are drawn 
in section [51 An auxiliary calculation of the elastic energy of a random dislocation wall 
is presented in the Appendix. 



2. Problem definition and simulations 



2.1. Equations of motion 

Consider a system of iV straight edge dislocations running parallel to the z axis of 
a Cartesian coordinate system. Let all dislocations have a common Burgers vector 
pointing along the x axis (a so-called single slip geometry), bj = Sj&e x , where Sj = ±1 is 
the sign of the zth dislocation. Assuming overdamped glide motion with a dislocation 
velocity v that is proportional to the local resolved shear stress, and zero dislocation 
mobility in the climb direction, the equation of motion of dislocation i piercing the xy 
plane at r*j = (xj, yi) can be written as 
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where \ denotes the dislocation glide mobility, G = p/[2ii(l — u)] where \i is the shear 
modulus and v is Poisson's ratio of the embedding isotropic crystal, r ind (r) denotes the 
resolved shear stress field induced by a positive dislocation located at the origin [15], 
and Text is a constant externally applied resolved shear stress. 

It is useful to introduce natural coordinates at this point which will be denoted by 
an apostrophe (') in the following. Measuring length in units of the average dislocation- 
dislocation distance p~ x l 2 (where p denotes the total dislocation density of dislocations 
including both signs and, in multiple slip geometries, including all slip systems), stress 
t in units of Gb^/p, and plastic strain 7 in units of b^fp leads to the relations 

x' = xy/p, t' = tpxGb 2 , t = t/ (Gby/p), 

v' = v/(J-p X Gb 2 ) } 7' = 7 /(Vp), Atf = AE/(Gb 2 ), 

where AE is the elastic energy difference between two states of the system (energy per 
unit dislocation length). In natural coordinates ([T]) takes the form 
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where (p denotes the angle between the x axis and r'. 
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2.2. Simulation details 

To study dislocation relaxation, a large number of discrete dislocation dynamics 
simulations have been performed. Equations of motion ([3]) were solved with the 4.5th 
order Runge-Kutta-Fehlberg method. Periodic boundary conditions were applied to a 
square simulation area with edges parallel to the slip planes, following the method used 
in [SJ. 

To avoid overly small timesteps during the final stages of approach of narrow dipoles 
(pairs of dislocations of opposite signs), a small number of extremely narrow dipoles 
were excluded from the solution of d3J) and forced to move as if they were isolated 
from the rest of the system. This is justified as the far- field stresses of the dislocations 
in a narrow dipole cancel, while their pair interaction diverges when the dislocation- 
dislocation distance approaches zero. As a consequence, the dynamics of narrow dipoles 
is effectively uncoupled from the rest of the dislocation system. We do not allow for 
annihilation of narrow dipoles, which is a process that is governed by the atomistics of 
the dislocation cores. This implies that we consider dislocation spacings to be large in 
comparison with the dipole annihilation distance, which is believed to be of the order 
of one nanometer. As typical dislocation densities in highly dislocated crystals are of 
the order of 10 14 m -2 , i.e. the average dislocation spacings are of the order of a hundred 
nanometers, this is not a severe restriction. 

A first set of simulations was started from random configurations of equal numbers 
of positive and negative dislocations. The number of simulated dislocations iV (which 
defines the system size L' s since N = Lf) varied between 16 and 128. It is well known 
that the flow stress of single-glide dislocation systems is around 0.1 in natural units 
[I]. To allow the dislocation systems to reach mechanical equilibrium at the end of the 
relaxation, we restricted the applied external stresses to levels below the flow stress, 
using stresses between and 0.088 natural units. 

As seen in figure HJ individual simulations showed strong avalanche-like activity 
during relaxation, as previously observed in [5]. To reveal scaling properties of the 
relaxation process, the evolution of global parameters such as stored energy, mean 
absolute dislocation velocity, mean square velocity, and mean strain rate was averaged 
over 3600 to 10 5 simulations starting from different random initial configurations. This 
ensemble averaging resulted in smooth ensemble averaged graphs as seen in figure [TJ 

In addition to the relaxation of 'neutral' arrangements of dislocations moving on a 
single slip system aligned parallel to the edge of a square simulation area, we considered 
two variants of this basic setting (these simulations were only performed at zero external 
stress): i) To elucidate the influence of different implementations of the boundary 
conditions, we performed single slip simulations with the simulation box oriented at an 
angle (p to the slip planes. This does not affect the short-range dislocation-dislocation 
interactions but modifies the stress field created by the periodic images, ii) To study 
the influence of a net Burgers vector on the relaxation process, we investigated the 
limiting case of fully polarised dislocation systems (dislocations of one sign only) of 
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Figure 1. Evolution of the average dislocation velocity for different ensemble sizes, 
sizes 16 < N < 128. 

To investigate the differences between single and multiple slip geometries, we also 
performed simulations in which dislocations of multiple slip systems were present. In 
these simulations we considered sets of equally populated slip systems, each containing 
edge dislocations with a zero net Burgers vector that were initially distributed at 
random. The methodology of the multiple slip simulations was identical to the single 
slip simulations described above, with the following differences: i) For each dislocation 
in the system, the complete elastic stress tensor was computed, again assuming periodic 
boundary conditions in a square simulation area and using the method of [16J. The forces 
acting on the dislocations were calculated from the stress tensor components using the 
Peach-Koehler formula [15] . ii) In addition to narrow dipoles, which were treated in the 
same way as described above, pairs of attracting dislocations on intersecting slip planes 
needed to be treated separately. Such dislocation pairs react with each other and form a 
reaction product ('dislocation lock') with Burgers vector equal to the net Burgers vector 
of the constituent dislocations. In real crystals, both the mobility of dislocation locks 
and the stress required for separating them into the constituent dislocations depend 
on their atomic core structure. For simplicity, we assumed all dislocation locks to be 
immobile and to possess infinite separation stress. New dislocations joining an existing 
lock were assumed to annihilate with the constituent dislocation of the opposite Burgers 
vector if such a dislocation was present. If a dislocation lock converted to an ordinary 
dislocation through such an annihilation event, it was assumed to become mobile again. 
To avoid overly small timesteps during the final stages of dislocation lock formation or 
reaction with a new dislocation (again due to diverging elastic interactions between the 
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constituent dislocations) the moving constituent dislocations were pinned when their 
relative distance decreased below a small predefined reaction radius which mimics the 
core extension of the dislocation lock. For simplicity, the stress field of a lock was 
calculated as the superposition of the stress fields of the pinned constituent dislocations. 
Although the positions of the constituent dislocations are scattered over a small region 
of the size of the reaction radius, their net stress field is a good approximation of the 
stress field of the lock at larger distances. Finally, we note that this stress field is of 
long range character since the net Burgers vector of a dislocation lock is not zero. 

3. Theory 

In this section a simple scaling argument is used to establish the asymptotic kinetics of 
a bimolecular combination reaction. We then adapt the fundamental ideas behind this 
argument to the relaxation of dislocation systems as described in the previous section. 
Based on the adapted model, predictions are made for the evolution of several physical 
quantities which can be directly obtained from the simulations described in section [2j 
These predictions are compared in detail to the numerical results in section HI 

3.1. Role of density fluctuations in bimolecular reaction kinetics 
Consider the direct combination reaction 



where c denotes the spatial average of the concentration field c(r). Note that the initial 
concentrations of molecules A and B are equal only in an average sense because of the 
thermally induced random positions of the individual molecules. These concentration 
fluctuations need to be taken into account when modelling the reaction kinetics. 
For situations where long-distance transport of particles occurs by free Brownian 
motion, a simple scaling argument which captures the essential physics was given by 
Ovchinnikov [T7] (which the reader is advised to consult for further details). Broadly 
speaking, this involves the subsequent dominance of two consecutive mechanisms 
directly corresponding to the different length scales inherent in the system: i) the 
reaction of those molecules that do not need long distance motion to find a reaction 
partner, followed by ii) the long range Brownian motion and subsequent reaction of 
excess reactants remaining as a result of initial concentration fluctuations. 

3.1.1. Stage 1: Reaction controlled kinetics. As discussed in [17] , the first stage of 
reaction (j3J) is controlled by a bimolecular reaction rate k which depends on the short- 
range interactions between the reaction partners. This stage can be discussed in the 
classical approximation 
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Solving dSJ) for initial conditions (J5J) yields 



2 \''((r,(r)-c[r) = rr;. 

(8) 



1 + fc^t (7) 
c A = c B = (kt)~\ t > (£;c A y 1 . 

3.1.2. Stage 2: Diffusion controlled kinetics. Because of thermal concentration 
fluctuations, not all of the molecules A and B can be consumed during stage 1. Spatial 
fluctuations of the reactant concentrations lead to local excess of molecules of one type. 
As these excess molecules cannot find a local partner, they need to migrate via long range 
Brownian motion, leading to a second kinetic stage controlled by long-range diffusion. 

To characterise the concentration fluctuations in question we observe that, on 
scales larger than the range of the 'contact interactions' which govern the first stage 
of the reaction kinetics, the positions of molecules are statistically independent. As a 
consequence, the initial numbers of molecules A or B in a sufficiently large volume V 
are Poisson distributed. This implies that the mean numbers of excess molecules fulfil 
the relations [17] 

Kr xC ess) 2 = ^(( CA (r)-cl) 2 ) = Vc A , 
(<Less) 2 = ^(( CB (r)-4) 2 ) = V4, 

where () denotes spatial averaging over a large statistically homogeneous system or, 
equivalently, ensemble averaging over a large number of statistically independent and 
equivalent realizations. 

To further discuss the reaction kinetics during stage 2, we confine ourselves to the 
limiting case of an infinite bimolecular reaction rate k — > oo. This choice does not affect 
the generality of the discussion, it only affects the moment of the crossover between 
stage 1 and stage 2 kinetics. We introduce the concentration difference z(r) as 

z(r) = c A (r) - c B (r). (9) 

In a hypothetical initial state before the reaction has been 'switched on' at t — 0, ca(p) 
and cb(V) are statistically independent and z(r) has the initial statistical properties [17] 

(z(r)) = (c A (r)-c B (r)) = 0, 

(z(rf) = <( CA (r) - c B (r)) 2 ) = 2<PjV = 2c° B /V. 
To proceed, we note that for k — > oo molecules A and B cannot coexist for t > 0. 
Therefore, z(r) gives a complete characterisation of the concentration map for t > 0: for 
z(r) > 0, ca(v) = z(r) and cb(t*) = and for z(r) < 0, ca(v) = and cb(V) = —z(r). 
Because there exists no other physical length scale in the system, the size of the regions 
characterised by z > and z < is determined by the diffusion length L(t) = \f~Dt (for 
simplicity, the same diffusion constant D is assumed for both kinds of reactants). 

Consider now the volume referring to the diffusion length at time t, V(t) = L(t) 3 = 
(Dt) 3 / 2 . One can suppose that for lengths larger than L{t) the fluctuations of z(r) are 
still not affected by diffusion; therefore, one can write that 

(\z(r)\) ~ (A/v{t)) l/2 = (c° A ) 1/2 
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In some regions, ca(p) = \z(r)\ and in others, cb(p) = hence, on average 

c A (t) = c B (£) = i(|z|> - ft) 1 * D-*'W, (12) 

meaning a slower kinetics than in stage 1 (CO). (For a more detailed derivation see |17j). 

Note that f|T2|) is only applicable within certain time limits [ti,^]- For instance, 
for k — > oo, t x is equal to the time the diffusion length L(t) needs to exceed the average 
intermolecular distance. £2 is determined by the time needed by L(t) to reach the system 
size L B , independent of the value of k. 

3.2. Role of density fluctuations in the relaxation of random dislocation configurations 

3.2.1. Analogy between chemical and dislocation systems. Now we consider the 
relaxation of initially random dislocation configurations with the same number of 
positive and negative dislocations (J2i s i — 0) following the equations of motion (j3J). 
Although not immediately evident, this relaxation process has strong phenomenological 
similarities to the kinetics of the chemical reaction described in section 13.11 

To elucidate the analogy, the relaxation process will be envisaged as a gradual 
screening of the long range elastic stress fields of individual dislocations through the 
formation of dislocation-dislocation correlations [9l [18]. We first envisage 'neutral' 
arrangements where dislocations of both signs are present in equal amounts. In a 
screened dislocation arrangement, the excess of one sign over the other has been 
eliminated on scales above a few dislocation spacings. Any dislocation arrangement 
where excess dislocations are completely eliminated can be envisaged as an assembly of 
dipoles where each dislocation has exactly one partner of opposite sign within a distance 
of the order of p _1//2 , and this picture will be used in the following argument. Hence, we 
envisage the relaxation process as the gradual formation of a large number of dislocation 
dipoles consisting of dislocations of opposite signs, i.e., as a bimolecular reaction process 
analogous to the chemical reaction (jlj). 

The principal difference between the two processes lies in the dynamics of individual 
particles: for the dislocation system, dislocation glide motion is driven by dislocation- 
dislocation interactions which scale like 1/r, whereas in case of the chemical reaction 
we are dealing with diffusive Brownian motion of the reactants. Despite this difference, 
we may again construct a two-stage model for the dislocation relaxation process: i) in 
stage 1 adjacent opposite sign dislocations form dipoles; ii) in stage 2 initial fluctuations 
in the excess dislocation density gradually die out from shorter towards longer length 
scales as excess dislocations which did not find a dipole partner in stage 1 undergo long 
range glide motion. The process terminates once the length scale on which fluctuations 
have been eliminated reaches the system size. 

3.2.2. Relaxation dynamics. To understand the time evolution of the system, we again 
consider the evolution of the typical length scale L(t) below which density fluctuations 
have already died out. To this end we consider that i) similar to the chemical case, areas 



Role of density fluctuations in the relaxation of random dislocation systems 



9 



of size L(t) 2 typically contain mobile excess dislocations (which have still not 'reacted' 
into dipoles) with only one or the other sign and that ii) dislocation dipoles give rise only 
to short range stress fields with a 1/r 2 decay. As the typical distance between opposite 
sign dislocations which try to find each other is proportional to L'(t'), the typical driving 
stress towards dipole formation scales as r' ~ 1/L'(t'). (Recall that variables with an 
apostrophe (') are measured in natural units.) The dislocation velocity is proportional 
to t', and therefore the characteristic time for eliminating excess dislocations on scale 
V scales as L'/t' oc {L') 2 . Hence we find that 
dL' 1 

or, equivalent ly L' ~ \p2±' . (13) 



dtf U 

Incidentally, this result is very similar to the evolution law of L(t) for Brownian motion, 
L = y/Dt. By supposing that the mobile dislocations inherit the initial concentration 
fluctuations we find that at time t', regions of size L'(t') 2 ~ 2t' contain about V excess 
dislocations of one or the other sign, and L' 2 dislocations in total (see the chemical 
model in section 13.1. 2D . Thus, the fraction of non-paired dislocations is estimated to 
decrease in time as 



-^excess L 1 



(14) 



N L' 2 yffi' 

Following the chemical model in section 13.1.21 it is straightforward to predict the 
time interval [t'j,^] during which the above argument is expected to hold. The start 
time t[ is characterised by L reaching the dislocation-dislocation distance p~ x l 2 (1 in 
natural units) and the process is finished when L reaches the system size L s . With f[T3T) 
this leads to 

1 2 2 2 2 2 t[ tt y 1 

in natural coordinates where A" denotes the total number of dislocations in the system. 



3.2.3. Scaling relations for energy and velocity. The fraction of 'non-paired' 
dislocations is not a convenient quantity for comparing the scaling model with 
dislocation dynamics simulations, as the definition of 'dislocation pairs' in a multipolar 
dislocation arrangement may be ambiguous. Instead, we consider the evolution of the 
nth moment (\v \ n ) of the dislocation velocity and of the excess elastic energy, both of 
which can be determined from the simulations in a straightforward manner. To obtain 
scaling estimates for the relaxation of these quantities, we assume that all dislocation 
dipoles are at rest and only the excess dislocations move. Furthermore, we assume that 
the motion of excess dislocations is not hindered by the dipoles already formed (this can 
be rationalised with the short range of the dipole stress field). The velocity of excess 
dislocations scales as \v'\ ~ l/L', leading to 

(wn - >< w\ n - ^ x = {LT [n+1) - (2tr (n+l)/2 , (i6) 

where ( 1T3"]) and f[T4"|) have been used. 
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The dynamics of dislocations is assumed to be overdamped. Hence, the work that 
is done by the internal stresses in driving the system is completely dissipated: the 
amount of dissipated energy exactly matches the reduction in elastic energy. The energy 
dissipated per unit time by a moving dislocation % scales like r^j (the Peach-Koehler 
force acting on the dislocation is proportional to the stress) , and consequently the time 
evolution of E can be expressed as 



Since the motion is overdamped, the dislocation velocity is proportional to the acting 
stress. In natural units and for t! e [^,^2], the ensemble averaged elastic energy thus 
evolves like 



where ( !T6|) was used to estimate the time evolution of the second velocity moment. 

3.2.4- Relaxation of arrangements of dislocations of the same sign. The scaling 
argument in the previous section is based on the formation of dipoles consisting of edge 
dislocations of opposite signs. At first glance, such an argument seems to be completely 
inapplicable to systems where only dislocations of the same sign are present. Dipole 
formation is clearly impossible in such systems. Instead, a most conspicuous feature in 
the relaxation of single-sign edge dislocation systems is the formation of walls containing 
many dislocations of the same sign that are aligned in the direction perpendicular to 
the slip plane [19]. Accordingly, theoretical arguments have focused on parameters 
characterising the 'condensation' of dislocations into walls |19j . However, even though 
dislocation wall formation is a most conspicuous feature, wall formation alone can not 
produce a screened dislocation arrangement. The authors of [19] evaluate the driving 
force for wall formation by assuming periodic spacing of dislocations along a wall and 
using classical results found e.g. in [15]. If we follow this line of reasoning and note that 
the energy of a dislocation in a wall decreases with decreasing dislocation spacing, the 
minimum-energy structure for the system at hand would be a single system-spanning 
wall. However, for an initially random dislocation system the y positions are independent 
random variables and it is not easy to see how a periodic arrangement could form in the 
absence of dislocation climb. We calculate the energy of a random wall in the Appendix 
and show that forming such a wall does not produce any energy reduction with respect 
to the initial random 2D arrangement. 

How then can an arrangement of dislocations of the same sign be screened? The 
answer was provided by Wilkens |8J who demonstrated that a screened dislocation 
arrangement can be constructed by eliminating dislocation density fluctuations above a 
certain scale. To this end, he proposed a construction where the crystal cross section 




i=i 



(17) 




(18) 



Role of density fluctuations in the relaxation of random dislocation systems 1 1 

is tiled into a grid of cells of size I, and the same number of dislocations is randomly 
distributed within each cell. This construction, which eliminates all density fluctuations 
on scales above the cell size, leads to an arrangement where the screening radius coincides 
with the cell size, S oc p\n(l/b). 

Taking the Wilkens construction as a well-screened reference state offers a surprising 
outlook on the relaxation of initially random systems of same-sign dislocations. With 
respect to this reference state, the initial random arrangement contains density 
fluctuations on all scales which may be either positive (5p = p(r) — (p) > 0, positive 
excess) or negative (Sp < 0, negative excess). To achieve screening, dislocations must 
migrate from regions of positive to regions of negative excess, and this process is governed 
by the long-range stress fields associated with the presence of (positive or negative) 
excess dislocations. In other words, the kinetics of the process follows from exactly 
the same scaling argument as used in the previous section: We are dealing with the 
stress-driven elimination of excess dislocation densities, with the only difference that 
the excess is now not of positive over negative dislocations, but of the local dislocation 
density over the average one. 

If the above argument is correct, the relaxation kinetics of same-sign dislocation 
systems should be characterised by a slow power-law stage which has the same 
characteristics as the relaxation of neutral dislocation systems as discussed in the 
previous sections. We demonstrate in the next section that this is indeed the case. 

4. Comparison with simulation results 

4-1- Single slip relaxation of dislocation arrangements with zero net Burgers vector 

4-1.1. Evolution of the elastic energy. The numerically determined evolution of the 
elastic energy is displayed in figure [2j The zero value of the energy was chosen to 
correspond to the final relaxed state of the system. As seen on the figure, the evolution of 
the elastic energy AE'/N per dislocation can be fitted satisfactorily with the prediction 
in (Tl8|) . AE/N = A/^/2F. Power-law relaxation occurs from times t[ ~ 0.3 onwards, 
in good agreement with the model prediction t[ ~ 1/2 in (fl5l) . The final equilibrium 
value of the elastic energy is not a priori known but was fitted to the data such as to 
achieve a maximum extension of the linear scaling regime. Unfortunately this precludes 
determination of the second critical time t' 2 from these results. The presence or absence 
of an external stress below the flow stress of the relaxed dislocation system seems to 
have negligible influence on the evolution of the elastic energy. 

Note that it is possible to collapse the curves for different system sizes L' s by 
normalising the graphs with ln(Lg), as done in figure El This observation is consistent 
with the fact that the elastic energy of the initial random dislocation system is of the 
order of Eq ~ NGb 2 ln(L s /6) while the energy of the final relaxed state is of the order of 
Eoo ~ NGb 2 In^^/pb) -1 ) [8]. For estimating the value of we used that the range 
of dislocation pair correlations in mechanical equilibrium is of the order of the mean 
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Figure 2. Evolution of the elastic energy per dislocation for different system sizes N 
and external stress values Tg Xt . 

dislocation-dislocation distance p~ 1 ^ 2 [9]. Therefore, the elastic energy difference per 
dislocation between the initial and final states is of the order of 

AE/N = {E - £oc)/iV ~ Gb 2 \n(L s ^p) = Gb 2 \n(L' s ). (19) 

Despite this relation connecting the initial and final states of the system, the numerical 
finding that the AE(t)/N curves for different system sizes can be collapsed on their 
entire course by normalising them with hi(Lg) is not trivial, as the agreement extends 
also to the relaxation kinetics and characteristic crossover time t[. 

4-1.2. Evolution of the mean square velocity. The numerically calculated evolution of 
the mean square velocity (v' 2 ) is displayed in figure[3]for zero applied stress and different 
system sizes. Due to the connection ( fT8l) between the mean square velocity and the 
elastic energy of the system, it is not surprising that similar statements apply here 
as for the evolution of the elastic energy. As seen in the figure, the model prediction 
(v' 2 ) ~ (2f)~ 3/2 in (HE]) fits the data well from t[ w 0.3. Due to the fact that (v' 2 ) 
is proportional to the time derivative of the elastic energy, its graphs are much noisier 
than those obtained for the energy, preventing again the detection of the supposed upper 
critical time t' 2 . As for the energy, size effects can be scaled out with a normalisation 
factor ln(Lg)" 1 which is a direct consequence of ([18]) and ( fT9l . A final analogy to the 
evolution of the elastic energy is that external stresses have only negligible influence on 
the evolution of the mean square velocity. For this reason, simulations with non-zero 
external stresses were omitted from figure El 
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Figure 3. Evolution of the mean square velocity (v 12 ) for different system sizes N at 
zero external stress. 

4-1.3. Evolution of the mean absolute velocity. In figure H] the evolution of the mean 
absolute velocity (\v'\) can be seen for different system sizes. Again, a power law 
time dependence can be observed from t[ ~ 0.3 although an exponent —0.86 gives 
a better fit than the theoretically predicted —1 expected according to equation ffTB"]) . 
One may argue that slowly moving dislocation dipoles play a bigger role in this case, 
as their small velocities contribute more strongly to (\v'\) than to (v' 2 ). Therefore, the 
gradually increasing number of slowly moving dislocations might be responsible for the 
reduced relaxation exponent. What makes this figure very interesting is the possibility 
to estimate values of t' 2 . It was found that t' 2 « 0.2iV gives a good approximation, in line 
with the model prediction t' 2 ~ N/2 in ( fT5l) . It was also observed that normalisation 
with (ln(Lg)) _1//2 collapses the graphs referring to different system sizes L' s in the region 
of small t' <C t' 2 . This is consistent with the relations for the elastic energy and the 
mean square velocity. Finally, as in case of the energy and the mean square velocity, 
the evolution of the mean absolute velocity is not changed by the presence of external 
stresses below the macroscopic flow stress. 

4-1.4- Evolution of the plastic strain rate. Another numerically measurable quantity 
is the plastic strain rate, defined as 

7 = E^ = &4- 2 E^ i' = N-^^[. (20) 

i S i i 

In the following the evolution of j 1 is studied for applied shear stresses r^ xt below 
the macroscopic flow stress for the present dislocation geometry. From the data of 
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Figure 4. Evolution of the mean absolute velocity (\v'\) for different system sizes N 
at zero external stress. 

Miguel and co-workers [U [6], this is estimated to be ~ 0.1 in natural units. As it was 
demonstrated in figure [2j external stresses in this range do not appreciably change the 
evolution of the elastic energy. We study the relaxation of the strain rate mainly in 
order to assess, by comparing with the work of Miguel et al., the relevance of different 
initial conditions on the creep behaviour of dislocation systems. 

Figure [5] shows the numerically determined evolution of the plastic strain rate j' for 
different system sizes and external stress values. As can be seen, the plastic strain rate 
scales roughly in proportion with the external stress. The relaxation does not follow any 
discernible power law but is roughly exponential. This is in marked contrast with the 
findings of Miguel et al. [H E] who for well relaxed initial configurations demonstrate 
an Andrade-type power-law decay, 7 oc t~ 2 ^ 3 . The discrepancy points to the crucial 
importance of initial conditions for relaxation processes in dislocation systems - a factor 
which is also borne out by the history dependence of creep relaxation processes that was 
demonstrated by Miguel et al. [6|. 

4-1.5. Role of boundary conditions (simulation box orientation). In this section we 
investigate the influence of different ways of implementing the periodic boundary 
conditions by tilting the angle between the edges of the simulation box and the trace of 
the slip planes. Simulations with different tilt angles <p are physically equivalent except 
for the spatial arrangement of the periodic images of each dislocation (see figure [6]). This 
arrangement affects the dislocation-dislocation interactions on scales comparable to the 
simulation box size. Also, the interaction energy of each dislocation with its periodic 
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Figure 5. Evolution of the plastic strain rate 7' normalised with the external stress 
Tg Xt for different system sizes N and external stress values r^ xt . 
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Figure 6. Arrangement of the first 8 periodic images of a given dislocation for 
<p = 0° (left) and ip = 45° (right). The latter configuration has a higher elastic 
energy as the nearest neighbours of each dislocation are in an energetically unfavourable 
configuration. 



images affects the initial elastic energy of the system, which is smallest for (p = and 
has a maximum for (p = 45°. 

The influence of simulation box orientation on the relaxation process is illustrated 
in figure The absolute values of the squared velocity (or equivalently the energy 
dissipation rate) are higher for ip = 45° than for ip = 0. However, both curves differ 
only by a constant factor (the ratio of the initial excess energies), while the dynamics 
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Figure 7. Evolution of the mean square velocity (v 12 ) in systems with a single slip 
geometry with the simulation box oriented at different angles ip to the slip planes. 

of the relaxation processes is otherwise identical. 

4-2. Relaxation of arrangements of dislocations of the same sign 

For systems of dislocations of the same sign we have numerically evaluated the evolution 
of the mean square velocity (v' 2 ) (or, equivalently, of the energy dissipation rate) and 
of the mean absolute velocity (\v'\). All calculations were performed at zero external 
stress since any applied stress would induce a sustained drift motion of the dislocation 
arrangement. 

Results are shown in figure M and figure [9] together with the fit functions obtained 
for neutral dislocation arrangements (see figure [3] and figure H]). It is evident that the 
relaxation of same-sign dislocation systems follows the same scaling laws that have 
been observed for systems containing equal numbers of dislocations of both signs. This 
provides strong support for our basic conjecture that relaxation is governed by the 
stress-driven elimination of excess dislocations in a process that progresses from small 
to large scales. The processes occurring on short scales, on the other hand, are evidently 
different for the two systems (dipole formation vs. formation of walls). This is reflected 
by the fact that the relaxation process in the single-sign dislocation systems shows an 
initial size dependence which is not present in neutral dislocation systems (see figure [3] 
and figure H] for comparison) . 
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Figure 8. Evolution of the mean square velocity (v 12 ) for different system sizes N for 
dislocations of the same sign. 




Figure 9. Evolution of the mean absolute velocity (\v'\) for different system sizes N 
for dislocations of the same sign. 
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Figure 10. Comparison of the evolution of the mean square velocity (v 12 ) in 
dislocation systems with single and multiple slip geometries. In the figure, p means 
the total density of dislocations on all slip systems. The last two curves were shifted 
downwards to help comparison. 



4-3. Dislocation relaxation on multiple slip systems 

Figure [TO] compares the relaxation of dislocation systems in single, double and triple 
slip. As seen in the figure, at long times the relaxation in multiple slip geometries 
accelerates in comparison with the relaxation in single slip. This is consistent with the 
idea that in multiple slip geometries relaxation processes proceed through the formation 
of dislocation locks and the annihilation of mobile dislocations at these locks. The long 
range stress fields associated with dislocation locks and the removal of dislocations 
accelerate the relaxation as seen in the figure. One of the main assumptions behind 
our scaling model, namely that the motion of mobile dislocations is governed mainly 
by their mutual interaction, no longer holds in multiple slip geometries. Therefore, the 
scaling model can not be applied to these situations. Indeed, as seen in figure [lOj no 
power law relaxation regime can be detected in the multiple slip simulations. 

5. Discussion and conclusions 

The present paper discusses the relaxation of initially random arrangements of straight, 
parallel edge dislocations. Following a phenomenological analogy with the kinetics of 
bimolecular reactions p2] the relaxation process can be divided into three consecutive 
stages. Stage 1 is characterised by rapid rearrangements of neighbouring dislocations, 
leading to the formation of dipoles, multipoles, and/or short wall segments. Stage 2 
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hosts the gradual extinction of initial fluctuations in the Burgers vector density on ever 
increasing length scales through the long range transport of excess dislocations. This 
stage gives rise to characteristic power-law relaxation dynamics. The relaxation process 
terminates in stage 3 when the characteristic fluctuation length reaches the system size. 

Our considerations focus on the power-law relaxation dynamics in stage 2. In 
section 3, we formulated a scaling theory for this process by making the analogy with 
a bimolecular reaction. To this end, we considered a highly simplified picture where 
pair interactions between positive and negative excess dislocations lead to long-range 
dislocation transport resulting in the formation of dislocation dipoles. This led to 
predictions for the evolution of the elastic energy and the first two moments of the 
dislocation velocity. These predictions were then compared to ensemble averages of 
discrete dislocation dynamics simulations, and convincing agreement was found for single 
slip geometries. For multiple slip geometries, however, the persistent long range stress 
fields of dislocation locks accelerate the relaxation process, to which the scaling model 
can no longer be applied. 

The actual dislocation processes in many- dislocation systems are much more 
complex than the simplified picture underlying our scaling arguments. Instead of 
long-range transport of excess dislocations and dipole formation, we see complex 
rearrangements resulting in dislocation dipoles, multipoles, and walls. In addition, it 
is well known that dislocations have a propensity to form large-scale heterogeneous 
patterns consisting of dislocation-rich and dislocation depleted regions. In the following 
we briefly discuss how these complex static and dynamic features fit into the idealised 
picture we used in the previous sections. 

5.1. Dipoles, multipoles and walls 

We have developed our scaling argument for dipole formation which can be considered 
a bimolecular reaction between positive and negative dislocations. However, actual 
dislocation arrangements are much more complex. It is therefore important to emphasise 
that the core of the argument is the elimination of large-scale fluctuations in the excess 
Burgers vector density, and not the resulting arrangement of nearby dislocations. For 
the long-time asymptotics of the relaxation process, which is governed by the elimination 
of fluctuations on larger and larger scales, the small-scale arrangement of dislocations 
is virtually irrelevant - at least as long as the local features (dipoles, multipoles, short 
walls and combinations of all these) do not give rise to long-range stresses. According to 
the present argument, a well-screened dislocation arrangement is one in which Burgers 
vector fluctuations have been eliminated on all scales. If we are dealing with a neutral 
dislocation system (equal numbers of positive and negative dislocations) this means that 
the net Burgers vector is zero in each small volume, i.e., for each dislocation we can find 
exactly one partner of opposite sign nearby. This motivates the dipolar picture even if 
the actual dislocation arrangements may be more complex. 

An alternative mechanism for creating well-screened dislocation arrangements is 
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the formation of system-spanning walls of dislocations of the same sign. Periodic 
arrangement of edge dislocations of one sign into a wall perpendicular to the glide 
plane removes the logarithmic divergence of the dislocation energy and introduces a 
screening length that is proportional to the dislocation spacing along the wall |15j . Walls 
are conspicuous both in simulations and in many experimentally observed dislocation 
microstructures. At first glance, wall formation mechanism seems to be completely at 
odds with the mechanism discussed in the present paper: Formation of walls of same- 
sign dislocations increases, rather than reduces, the Burgers vector density fluctuations. 
However, a closer investigation reveals that wall formation by itself is not a screening 
mechanism at all. Forming a wall of randomly spaced dislocations does not reduce 
the energy in comparison with a random 2D dislocation arrangement (see Appendix). 
Instead, the screening effect is contingent on the equal spacing of dislocations, i.e. on 
suppressing fluctuations of the Burgers vector density along the wall direction. If we 
start from a random dislocation arrangement this is not easy to obtain: Either the 
dislocations must have climb degrees of freedom (which we do not consider in the 
present study), or dislocation motions that lead to the formation of multiple walls must 
be correlated over large distances such as to ensure that each wall collects only those 
dislocations that fit into an evenly spaced pattern. In the latter case, we are again 
dealing with the suppression of Burgers vector density fluctuations on all scales above 
the wall spacing, and the long time asymptotics of this process is expected to obey our 
scaling theory. This is confirmed by the simulations, which however also demonstrate 
that the short-time behaviour is different for neutral dislocation systems where the 
local arrangement of dislocations is characterised by dipolar and multipolar patterns, 
and for single-sign dislocation systems where the local arrangement of dislocations is 
characterised by walls (compare figure [3] and figure [8]). 

5.2. Long-range transport of excess dislocations 

Our scaling argument considers the stress-driven long-range transport of excess 
dislocations. The picture underlying the argument is schematically shown in figure [11] 
(top): A positive excess dislocation at A is attracted by a negative excess dislocation at 
B and the two recombine by long-range motion which is not affected by the stress field 
of the dislocations in between A and B. The figure also indicates that this idealisation 
may not be feasible when we are dealing with multipolar arrangements rather than 
isolated narrow dipoles: In that case the mutual interaction of the excess dislocation 
may be much weaker than their interaction with other dislocations 'on the way'. As a 
consequence, recombination is much more likely to occur by a collective rearrangement 
as shown in figure [11] (bottom). 

How does this affect our scaling argument? The total driving force for the process 
is the same in both cases. However, a collective rearrangement on scale L is likely to 
involve N = L^fp dislocations [N = 3 in figure [TT1 (bottom)]. Hence the driving force 
per dislocation is reduced by a factor of the order of N. However, the same is true for 
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Figure 11. Schematic illustration of excess dislocation transport in a multipolar 
dislocation arrangement: a positive excess dislocation at A 'recombines' with a negative 
excess dislocation at B (top) by direct motion from A to B or (bottom) by collective 
rearrangement involving multiple dislocations. 



the characteristic distance that has to be covered by each dislocation (L in the case of 
direct transport, 1/a/p in the case of collective rearrangement). As we assume that the 
dislocation velocity is proportional to the driving force, it follows that the characteristic 
time scale for eliminating the excess dislocation is the same in both cases, and our 
scaling argument remains valid. 



5.3. Large-scale dislocation patterning 

It is a well known phenomenon that dislocation microstructures forming during 
plastic deformation form heterogeneous patterns consisting of regions of high and 
low dislocation density, with characteristic lengths that are large in comparison with 
the dislocation spacing. On a conceptual level, possible mechanisms underlying this 
patterning were discussed by Nabarro [20] who pointed out that it may be energetically 
favourable to 'segregate' the dislocation microstructure into areas of high and low 
density: If we are dealing with a well-screened dislocation arrangement, the energy 
density scales like £ oc pln(l/(6 v /p)). In this case it can be easily shown that it is 
energetically favourable to increase p in some regions and decrease it proportionally 
in others. Such 'energetically driven' dislocation patterning could be a reason for the 
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formation of dislocation-dense and dislocation-depleted regions that is observed in many 
experiments. While this mechanism is not covered by the present model, it is not at 
variance with our considerations: In a neutral dislocation arrangement, the formation of 
dislocation-dense and dislocation-depleted regions might occur without disturbing the 
Burgers vector balance. We note, however, that in our simulations large-scale dislocation 
patterning is not observed - either because it does indeed not occur in single slip, or 
because the dislocation numbers in our simulations might be too small. 

5.4- Conclusion 

In conclusion we discuss the relevance of the processes discussed in the present paper 
for real-world systems. The relaxation of a random dislocation system has no direct 
counterpart in real deformation experiments, since it is impossible to 'prepare' such a 
random system in the first place. Our analysis of the screening of same-sign dislocation 
systems is, however, of general importance for understanding real dislocation patterns 
since it demonstrates that wall formation, though a conspicuous feature, can by itself 
not account for screening. This observation points to the importance of investigating 
long-range correlations between dislocation positions both within the walls and across 
different walls, and offers ample scope for future investigations. 

The investigated processes are of significant importance for discrete dislocation 
dynamics simulations of plasticity as the slow nature of the relaxation makes it difficult 
to obtain well-defined and energetically stable initial configurations. Our comparison 
of strain-rate relaxation experiments with the results of Miguel et al. demonstrates 
that the collective behaviour of dislocation systems may depend significantly on initial 
conditions. The analysis of this dependence is still in its infancy, yet understanding it 
is indispensable for carrying out dislocation plasticity simulations in a controlled and 
well-defined manner. 
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Appendix: Calculation of the elastic energy of a random dislocation wall 

We consider a wall of infinite height running along the plane x = in an isotropic 
material. Edge dislocations of Burgers vector b = be x are distributed randomly along 
the wall with average linear density 1/h. The geometry corresponds to a plane-strain 
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situation, hence the elastic energy density can be written as 

The total energy of the system is obtained by integrating (IA.1I) over the system volume, 

E = J v £d 2 r = |((T W + o-yy) 2 d 2 r, (A.2) 

where we have used that, for an infinite system, the second and third terms on the 
right-hand side do not contribute to the total energy. This can be shown as follows: For 
plane-strain deformation, the stresses can be written as derivatives of the Airy stress 
function x, ° xx = d 2 x, v yv = <%X anci &xy = ~d x d y x- Hence, 

J Wly ~ o- xx a yy }d 2 r = J v \d x d y xd x d y x - d 2 x xd 2 y x]& 2 r. (A.3) 

Partially integrating the second term in the integral on the right-hand side with respect 
to x and y shows that this integral contributes only surface terms to the total energy. 
These terms are negligible in the infinite- system limit. 

The ensemble-averaged stress at any point is given by summing over the stress fields 
of the individual dislocations in the wall and averaging over the different realizations of 
the random dislocation positions: 

Mr)) = (r)>, (A-4) 

n 

where aff is the ij component of the stress created at r by the nth dislocation. The 
elastic energy of the system depends on the averages of products (cfij(r)aki(r)) where 
(ij,kl) G [(xx, xx); (xx, yy); (yy, yy)]. In evaluating these averages we use that the y 
coordinates of the individual dislocations are independent random variables: 

{<Tij(r)<T M (r)) = (<r -(r))(<T W (r)) 

+ E [<4V)4V)> - (4 ) (r))(4? ) (r))] . (A.5) 

n 

We now make the following observations: 

• The average stresses {o~ij(r)) and their products (aij(r)aki(r)) depend on the x 
coordinate only. 

• The average single-dislocation stresses (c^) and (c^) become zero in the limit 
a — > oo, since these stresses are antisymmetric functions of the y coordinate. The 
same is true for the average total stresses (a xx ) and (a yy ). 

With these observations and using ( 1A.5I) we can write the system energy as 

NGb 2 n(l-u) In f L 



m ■ (A - 6) 

where the second step follows by interchanging the averaging and the integration. L 
is the system size (tending to infinity) which, in the absence of any other screening 
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mechanism, delimits the divergence of the dislocation self-energy. For periodic boundary 
conditions, as used in our simulations, L must be understood as the size of the 
periodic simulation box which in this case defines the screening length for an otherwise 
uncorrelated dislocation arrangement. iV = L/h is the total number of dislocations in 
the system. 

It follows from flA.61) that the energy per dislocation is equal to the energy of a 
single unscreened dislocation and, hence, equals the energy in a completely random 2D 
arrangement. In other words, the arrangement of dislocations of the same sign in a 
random wall does (with the possible exception of surface terms that are negligible in 
the infinite system limit) not produce any reduction of the total energy. Accordingly, 
the thermodynamic driving force towards forming such a wall is zero. 
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